{ "cells": [ { "cell_type": "markdown", "id": "7030540f", "metadata": {}, "source": [ "# Multiple Grid Cells in MUSICA\n", "In MUSICA, the State object that is present within the Solver object has an attribute called number_of_grid_cells.\n", "This attribute dictates the number of independent sets of well-mixed air masses whose chemical system will be solved by the same numerical solver.\n", "This tutorial will go over a simple example of solving a multi-grid-cell chemical system in MUSICA." ] }, { "cell_type": "markdown", "id": "d7d9a125", "metadata": {}, "source": [ "## MUSICA: Before Getting Started\n", "\n", "It is heavily recommended to go through the tutorials in the MusicBox repository first before going through the MUSICA ones since the former explains the workflow used in the latter.\n", "However, if you are only interested in MUSICA, here's how you can set up a virtual environment for it:\n", "\n", "```\n", "conda create --name musica python=3.9\n", "conda activate musica\n", "pip install musica\n", "conda install ipykernel scikit-learn seaborn scipy dask\n", "```" ] }, { "cell_type": "markdown", "id": "008e870b", "metadata": {}, "source": [ "## 1. Importing Libraries\n", "Below is a list of the required libraries for this tutorial:" ] }, { "cell_type": "code", "execution_count": 1, "id": "7c921c61", "metadata": {}, "outputs": [], "source": [ "import musica\n", "import musica.mechanism_configuration as mc\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import numpy as np\n", "pd.set_option('display.float_format', str) # This is done to make the arrays more readable\n", "np.set_printoptions(suppress=True) # This is done to make the arrays more readable" ] }, { "cell_type": "markdown", "id": "01472a4f", "metadata": {}, "source": [ "## 2. Defining a System\n", "\n", "This code snippet is a MUSICA version of setting up a system, which has an identical workflow to MusicBox.\n", "For an explanation of this code in MusicBox, please refer to the Basic Workflow tutorial in the MusicBox repository." ] }, { "cell_type": "code", "execution_count": 2, "id": "6b1084ee", "metadata": {}, "outputs": [], "source": [ "A = mc.Species(name=\"A\") # Create each of the species with their respective names\n", "B = mc.Species(name=\"B\")\n", "C = mc.Species(name=\"C\")\n", "species = [A, B, C] # Bundle the species into a list\n", "gas = mc.Phase(name=\"gas\", species=species) # Create a gas phase object containing the species\n", "\n", "r1 = mc.Arrhenius( # Create the reactions with their name, constants, reactants, products, and phase\n", " name=\"A_to_B\",\n", " A=4.0e-3, # Pre-exponential factor\n", " C=50, # Activation energy (units assumed to be K)\n", " reactants=[A],\n", " products=[B],\n", " gas_phase=gas\n", ")\n", "\n", "r2 = mc.Arrhenius(\n", " name=\"B_to_C\",\n", " A=4.0e-3,\n", " C=50, \n", " reactants=[B],\n", " products=[C],\n", " gas_phase=gas\n", ")\n", "\n", "mechanism = mc.Mechanism( # Define the mechanism which contains a name, the species, the phases, and reactions\n", " name=\"musica_micm_example\",\n", " species=species,\n", " phases=[gas],\n", " reactions=[r1, r2]\n", ")" ] }, { "cell_type": "markdown", "id": "7861e510", "metadata": {}, "source": [ "## 3. Creating the Solver\n", "Something more unique to MUSICA is that you have to manually define the solver.\n", "A solver integrates the chemical reactions that determine how atmospheric chemistry proceeds over time.\n", "There are a handful of solvers available, but Rosenbrock Standard Order will be used here.\n", "For more information on the types of solvers available, go [here](https://ncar.github.io/micm/user_guide/solver_configurations.html)." ] }, { "cell_type": "code", "execution_count": 3, "id": "18bb7d49", "metadata": {}, "outputs": [], "source": [ "solver = musica.MICM(mechanism = mechanism, solver_type = musica.SolverType.rosenbrock_standard_order)" ] }, { "cell_type": "markdown", "id": "2199a26f", "metadata": {}, "source": [ "## 4. Creating the State\n", "For MUSICA, the state must be created manually as well, with the number of grid cells being passed into the create_state() function.\n", "The state represents everything pertinent to solve chemistry. By solving, we meaning determining the concentrations at the next time step.\n", "Feel free to change the num_grid_cells value to experiment yourself." ] }, { "cell_type": "code", "execution_count": 4, "id": "d91a4be8", "metadata": {}, "outputs": [], "source": [ "num_grid_cells = 2\n", "state = solver.create_state(num_grid_cells)" ] }, { "cell_type": "markdown", "id": "54f3e880", "metadata": {}, "source": [ "## 5. Populating the Grid Cells\n", "5 dimensions will be used to populate the data for each of the air masses:\n", "* temperature (Kelvin),\n", "* pressure (Pascals), and\n", "* the concentrations of each of the species (mol/m3).\n", "\n", "The NumPy array has the reshape() function called on it so that it can be split up in the next step.\n", "The two arguments specify the number of rows and columns of the array, where -1 means any number of rows.\n", "Do note that the ordering inside the array matters and cannot be changed." ] }, { "cell_type": "code", "execution_count": 5, "id": "1dadda9c", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 300. , 101253.3, 5. , 5. , 5. ],\n", " [ 100. , 11253.3, 20. , 3. , 7. ]])" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "box_model_values = np.array([[300, 101253.3, 5, 5, 5], [100, 11253.3, 20, 3, 7]])\n", "box_model_values = box_model_values.reshape(-1, 5)\n", "display(box_model_values)" ] }, { "cell_type": "markdown", "id": "08d6e59c", "metadata": {}, "source": [ "## 6. Splitting up the Array Output\n", "Next, the values from the box_model_values array are taken and populated into variables so that they can be passed into the solver's state.\n", "The sample is organized into 5 columns that represent the 5 variables.\n", "The sample's rows each represent a distinct grid cell.\n", "These columns are populated into their respective variables and then passed into the solver's state.\n", "Do note that set_conditions() and set_concentrations() must be called with their respective arguments for the solver to successfully run.\n", "The concentrations must go through the extra step of being bundled into a dictionary since MUSICA explicitly requires a dictionary argument in the set_concentrations function.\n", "Lastly, an empty array is initialized to represent the solved concentration array at every time step, as well as the time step length, simulation length, and the current time step (all in seconds)." ] }, { "cell_type": "code", "execution_count": 6, "id": "e0a74a68", "metadata": {}, "outputs": [], "source": [ "temperatures = box_model_values[:, 0]\n", "pressures = box_model_values[:, 1]\n", "concentrations = {\n", " \"A\": [],\n", " \"B\": [],\n", " \"C\": []\n", "}\n", "concentrations[\"A\"] = box_model_values[:, 2]\n", "concentrations[\"B\"] = box_model_values[:, 3]\n", "concentrations[\"C\"] = box_model_values[:, 4]\n", "\n", "state.set_conditions(temperatures, pressures)\n", "state.set_concentrations(concentrations)\n", "concentrations_solved = []\n", "time_step_length = 1\n", "sim_length = 60\n", "curr_time = 0" ] }, { "cell_type": "markdown", "id": "1238f02b", "metadata": {}, "source": [ "## 7. Running the Solver\n", "This code solves the system at every specified time step and the solved concentrations are appended to the array.\n", "The first step will always be the initial conditions since at time = 0 seconds the reaction has not begun." ] }, { "cell_type": "code", "execution_count": 7, "id": "76b21308", "metadata": {}, "outputs": [], "source": [ "while curr_time <= sim_length:\n", " solver.solve(state, curr_time)\n", " concentrations_solved.append(state.get_concentrations())\n", " curr_time += time_step_length" ] }, { "cell_type": "markdown", "id": "5d2185c0", "metadata": {}, "source": [ "## 8. Preparing the Results (Advanced; Optional Read)\n", "Here, a new array is made that grabs only the first value (grid cell) for each key (species) at every time step.\n", "The time step index is divided by the time_step_length to account for time step lengths that are greater than one for proper array indexing.\n", "That new array is then passed into a Pandas DataFrame with the concentration columns renamed.\n", "Next, a time column is created from a range that represents the elapsed time at each time step.\n", "In this simulation, the temperature, pressure, and air density are all constant, so numpy's repeat() function is used to repeat their respective values for every time step.\n", "Once all the attributes are added to the DataFrame, their order is changed to follow a more logical flow.\n", "Due to the complexity of this code cell, it has been bundled into a function that takes in an argument for the grid cell you wish to convert to a DataFrame.\n", "The function is then called twice, one for each grid cell index." ] }, { "cell_type": "code", "execution_count": 8, "id": "65de1f01", "metadata": {}, "outputs": [], "source": [ "def convert_results_single_cell(cell_index):\n", " concentrations_solved_pd = []\n", " for i in range(0, sim_length + 1, time_step_length):\n", " concentrations_solved_pd.append({species: concentration[cell_index] for species, concentration in concentrations_solved[int(i/time_step_length)].items()})\n", " df = pd.DataFrame(concentrations_solved_pd)\n", " df = df.rename(columns = {'A' : 'CONC.A.mol m-3', 'B' : 'CONC.B.mol m-3', 'C' : 'CONC.C.mol m-3'})\n", " df['time.s'] = list(map(float, range(0, sim_length + 1, time_step_length)))\n", " df['ENV.temperature.K'] = np.repeat(temperatures[cell_index], sim_length/time_step_length + 1.0)\n", " df['ENV.pressure.Pa'] = np.repeat(pressures[cell_index], sim_length/time_step_length + 1.0)\n", " df['ENV.air number density.mol m-3'] = np.repeat(state.get_conditions()['air_density'][cell_index], sim_length/time_step_length + 1.0)\n", " df = df[['time.s', 'ENV.temperature.K', 'ENV.pressure.Pa', 'ENV.air number density.mol m-3', 'CONC.A.mol m-3', 'CONC.B.mol m-3', 'CONC.C.mol m-3']]\n", " return concentrations_solved_pd, df" ] }, { "cell_type": "code", "execution_count": 9, "id": "e9dda1ed", "metadata": {}, "outputs": [], "source": [ "concentrations_solved_pd_0, df_0 = convert_results_single_cell(0)\n", "concentrations_solved_pd_1, df_1 = convert_results_single_cell(1)" ] }, { "cell_type": "markdown", "id": "ada2398b", "metadata": {}, "source": [ "## 9. Viewing the Results\n", "With the DataFrames being fully prepared now, they are displayed and plotted to show the evolution of both of the systems over time." ] }, { "cell_type": "code", "execution_count": 10, "id": "d49f1b08", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", " | time.s | \n", "ENV.temperature.K | \n", "ENV.pressure.Pa | \n", "ENV.air number density.mol m-3 | \n", "CONC.A.mol m-3 | \n", "CONC.B.mol m-3 | \n", "CONC.C.mol m-3 | \n", "
---|---|---|---|---|---|---|---|
0 | \n", "0.0 | \n", "300.0 | \n", "101253.3 | \n", "40.59324282282551 | \n", "5.0 | \n", "5.0 | \n", "5.0 | \n", "
1 | \n", "1.0 | \n", "300.0 | \n", "101253.3 | \n", "40.59324282282551 | \n", "4.976428528398802 | \n", "4.999944351093875 | \n", "5.0236271205073235 | \n", "
2 | \n", "2.0 | \n", "300.0 | \n", "101253.3 | \n", "40.59324282282551 | \n", "4.929618429433882 | \n", "4.999502304184354 | \n", "5.070879266381759 | \n", "
3 | \n", "3.0 | \n", "300.0 | \n", "101253.3 | \n", "40.59324282282551 | \n", "4.860227571898497 | \n", "4.998027906909858 | \n", "5.141744521191641 | \n", "
4 | \n", "4.0 | \n", "300.0 | \n", "101253.3 | \n", "40.59324282282551 | \n", "4.769223461910271 | \n", "4.994590343685047 | \n", "5.236186194404681 | \n", "
... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
56 | \n", "56.0 | \n", "300.0 | \n", "101253.3 | \n", "40.59324282282551 | \n", "0.002651107160467822 | \n", "0.022649502677897896 | \n", "14.974699390161621 | \n", "
57 | \n", "57.0 | \n", "300.0 | \n", "101253.3 | \n", "40.59324282282551 | \n", "0.0020250213465128905 | \n", "0.01784640325737679 | \n", "14.980128575396096 | \n", "
58 | \n", "58.0 | \n", "300.0 | \n", "101253.3 | \n", "40.59324282282551 | \n", "0.0015394935543717359 | \n", "0.01398971727457048 | \n", "14.984470789171041 | \n", "
59 | \n", "59.0 | \n", "300.0 | \n", "101253.3 | \n", "40.59324282282551 | \n", "0.0011648552812954611 | \n", "0.010910311593626484 | \n", "14.987924833125062 | \n", "
60 | \n", "60.0 | \n", "300.0 | \n", "101253.3 | \n", "40.59324282282551 | \n", "0.0008772265708207769 | \n", "0.008465235540415044 | \n", "14.990657537888747 | \n", "
61 rows × 7 columns
\n", "\n", " | time.s | \n", "ENV.temperature.K | \n", "ENV.pressure.Pa | \n", "ENV.air number density.mol m-3 | \n", "CONC.A.mol m-3 | \n", "CONC.B.mol m-3 | \n", "CONC.C.mol m-3 | \n", "
---|---|---|---|---|---|---|---|
0 | \n", "0.0 | \n", "100.0 | \n", "11253.3 | \n", "13.53460893002309 | \n", "20.0 | \n", "3.0 | \n", "7.0 | \n", "
1 | \n", "1.0 | \n", "100.0 | \n", "11253.3 | \n", "13.53460893002309 | \n", "19.868536268722355 | \n", "3.1113111545426153 | \n", "7.020152576735027 | \n", "
2 | \n", "2.0 | \n", "100.0 | \n", "11253.3 | \n", "13.53460893002309 | \n", "19.608195525999147 | \n", "3.329170747441501 | \n", "7.062633726559351 | \n", "
3 | \n", "3.0 | \n", "100.0 | \n", "11253.3 | \n", "13.53460893002309 | \n", "19.22406658832665 | \n", "3.6442931108015326 | \n", "7.131640300871819 | \n", "
4 | \n", "4.0 | \n", "100.0 | \n", "11253.3 | \n", "13.53460893002309 | \n", "18.723574920030497 | \n", "4.043334633782276 | \n", "7.233090446187235 | \n", "
... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
56 | \n", "56.0 | \n", "100.0 | \n", "11253.3 | \n", "13.53460893002309 | \n", "0.0005360881626445901 | \n", "0.005726115410867138 | \n", "29.993737796426494 | \n", "
57 | \n", "57.0 | \n", "100.0 | \n", "11253.3 | \n", "13.53460893002309 | \n", "0.00036805009976445557 | \n", "0.0040698479741160445 | \n", "29.995562101926122 | \n", "
58 | \n", "58.0 | \n", "100.0 | \n", "11253.3 | \n", "13.53460893002309 | \n", "0.00025101920969299854 | \n", "0.002871932602189201 | \n", "29.996877048188125 | \n", "
59 | \n", "59.0 | \n", "100.0 | \n", "11253.3 | \n", "13.53460893002309 | \n", "0.00017007318795390593 | \n", "0.0020121301960677983 | \n", "29.997817796615987 | \n", "
60 | \n", "60.0 | \n", "100.0 | \n", "11253.3 | \n", "13.53460893002309 | \n", "0.0001144703965387078 | \n", "0.0013996875381660917 | \n", "29.998485842065303 | \n", "
61 rows × 7 columns
\n", "